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The effect of three-body interatomic contributions in the equation of state of He are investigated. 

^Nj , A recent two-body potential together with the Cohen and Murrell (Chem. Phys. Lett. 260, 371 

j_l ' (1996)) three-body potential are applied to describe bulk helium. The triple-dipole dispersion and 

^ l' exchange energies are evaluated subjected only to statistical uncertainties. An extension of the 

, diffusion Monte Carlo method is applied in order to compute very small energies differences. The 

■ results show how the three-body contributions affects the ground-state energy, the equilibrium, 

' melting and freezing densities. 
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I. INTRODUCTION 



O ' The unique properties of the helium systems at low temperature have attracted a continuous experimental and 
' theoretical interest in the investigation of their ground state potential energyi In the past, the construction of the 
• ' best potentials used semi-empirical methods where some parameters were obtained by fits to experimental data. One 
^ I of them, the so called HFDHE2 potential of Aziz and collaborators has allowed the understand of many properties 
• ^ . of helium in the condensed phases^*^ Despite of small inconsistencies in this potential, it was used during a long time 
in these studies. It was attractive to use an effective pair-wise additive potential and avoid considering high-order 
(-H , interactions among the atoms. 

O i' In the last decade, after a bound helium dimer was observed;^ great efforts were applied to develop ab initio 
methods in the description of He-He potentials This approach was very successful; interacting energies calculated 
^ • using infinite order symmetry adapted perturbation theory. Green's function Monte Carlo results and accurate dispcr- 
' sion coefficients were fitted to a Tang-Tonnies model^ and produced^ to date the best characterization of the helium 
. potential energy>i£ 

^-H ' The use of very accurate two-body potential energies, like those offered by the ab initio potentials, in the investi- 
^-H . gation of bulk helium, uncovered what was known for a long time; the correct description of many of its properties 
\l ' requires more general many-body potential. Among the many recent works where this situation was observed, we cite 
O ; Refs. irilll^ . 

• The first attempt to incorporate long range three-body interactions were made by Axilrod- Teller— and Muto— 
They used third order perturbation theory to calculate the triple-dipole dispersion energy for spherical symmetric 
^ ■ atoms. Effects of three-body exchange of electrons in trimers of helium started 10 years latter by Roseni^ using a 
valance bond approach. Since the work of Jansen and collaboratorsiii^ it is conjectured that three-body exchange 
energy is needed to understand the energy difference between the fee and hep crystalline structure. Many developments 
(-H ' occurred in the investigation of nonadditive effects are reviewed in Ref.[l9j. 

It is known that two-body interactions favor the hep over the fee structure. The hep crystalline lattice is about 
1!^ ■ 0,01% more stable than the corresponding fee structure— The two-body interaction potentials of the rare gas in 
general are very well known and the residuals errors associated to them can not be responsible for discrepancies 
in the calculated ground state energy and in the favored crystalline structure. At low temperatures all rare gas 
^ ] solids, but helium, crystallizes in a fee lattice. To overcome these discrepancies higher order interactions need to 
be considered. The inclusion of triple-dipole interaction improved the agreement between the experimental and 
theoretical values of the ground state energy, but still left the hep as the more estable structure. The inclusion of 
higher order terms in the dispersion energy like dipole-octopole and quadrupole-quadrupole terms did not modify 
the theoretical results significantly. So far, the proposed higher-order interactions greater than three-body have given 
no significant contributions to the interatomic potential. In this context it is important to better understand the 
three-body exchange interaction and have at the same time reliable ways to compute its associate energy. 

As already mentioned, for the helium systems the largest and most well know part of the three-body interatomic 
potential is the triple-dipole term of Axilrod-Teller-Muto. The exchange contribution is less known. It is much more 
weak than the dispersion energy and might be of importance in the crystallization process. Despite of its importance, 
still now, competing calculations might differ by an order of magnitudeSi and different potentials forms are used for 
fitting theoretical results obtained in calculations of the exchange energy. 
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Quantum Monte Carlo methods, where quantities of interest can be computed subject only to statistical uncer- 
tainties, can be very useful in the analysis and understanding of the different contributions to the potential energy. 
However straightforward use of these methods to compute small energies differences might not be possible. Results 
of independent runs and their associate statistical uncertainties might render a simple energy difference meaningless. 
Perturbation methods used together with quantum Monte Carlo methods still needs extrapolation^ that introduces 
further uncertainties. 

To be able to compute energies subject only to statistical uncertainties (avoiding extrapolation) is not only a 
matter of principle but a necessity in the present case, where we have a delicate balance between energies. A better 
understanding of the individual contributions of the three-body interactions are important by themself and moreover 
can increase the physical content of the analytical functional forms used to fit their contributions. 

In this work, we have applied a recent two-body potential^ together with the three-body potential of Cohen and 
MurrellS^ to describe properties of bulk ^He. By developing an extension of the diffusion Monte Carlo (DMC) method 
we are able to compute and analyze the individual contributions of the Coulomb and exchange terms of the three-body 
interactions. The effect of these contributions are considered in the equation of state and quantitative results show 
how they affect the ground-state energy, the equilibrium, melting and freezing densities. 

The paper is organized as follows, in the next section we present the Hamiltonian together with the interacting 
potentials used in this work. In Section IIIII we briefly describe one of the standard implementation of the Diffusion 
Monte Carlo algorithm, then we present an extension of this algorithm. It allows the use of a single set of walkers and 
reweighting to compute properties of a system of helium atoms described by different interacting potentials. Section 
[^contains details of our simulations; the results are presented in Section A discussion in Section IVll concludes 
the work. 



II. THE MODEL 



The Hamiltonian we use to describe the system of helium atoms is given by 



" = -^^n + V{R), (2.1) 

where R = {ri, r2, . . . , tn} stands for the N coordinates of the helium atoms and V{R) is the interatomic potential. In 
this work three sets of calculations were performed. In the first one the interatomic potential employed is an additive 
pair- wise potential V2{R) as proposed by Aziz and co-worksi^ In a second set, we considered the V2D potential 

V2D{R)^V2iR) + VD{R), (2.2) 
obtained by adding to V2 a damped Axilrod-Teller-Muttoi^*i^*2^ triple-dipole term (ddd = D) 



rr ry(3) 1 1 ^ ^ + ^ COS 7l COS 72 COS 73 

Vd = Z^-'''(lll) -F(ri2,ri3,r23), (2.3) 

[ri2ri3r23r 

where Z('^^(lll) is a constant, the ji are the internal angles of the triangle (formed by the three particles) and the 
rij the lengths of its sides. The damping F is given by the product 



-F'(?'i2, ri3, r23) = /(ri2)/(ri3)/(r23) (2.4) 

that depends on 

^^"^^^^"U otherwise, ^^'^^ 

where / and k are parameters. The damping of the dispersion energies at distances where charge overlap is significant 
is need for a reasonable description of the short range forces. The value used for ^^^^(111) in Eq. I|2.3|l is 0.324 K, as 
obtained by double perturbation theory^ 

Finally, the most complete interatomic potential we have considered. 



V2Dj{R) = V2{R) + Vn{R) + Vj{R), 



(2.6) 
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includes contributions from the exchange potential Vj as well. The two last terms of Eq. (|2.6|l form the three-body 
potential proposed by Cohen and MurrelliS^ The potential Vj is expressed through symmetry adapted coordinates Qi 
(linear combinations of the three distances r^j) 



Its functional form is given by: 



Qi = "^(''is + ri3 + r23), 

Q2 = -^(rig - ras), 

Q3 = ^(2ri2 - ri3 - r23). 
V6 



(2.7) 
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III. THE DIFFUSION MONTE CARLO METHOD 



A. The Standard Algorithm 



In almost all practical implementations of the diffusion Monte Carlo method2S*S& we compute quantities of interest 
by sampling the probability distribution fo{R) — 4'g{R)'<Po{R) that depends on V'o(-R), the true ground state wave 
function of the system and on a given guiding function, ^jq{R). 

A careful version of this method is presented by Umrigar, Nightingale and RungejSl A slight simpler implementation 
could be described as follow. It is convenient to start the calculation with a set of configurations draw from iV'Gp, 
obtained through the Metropolis algorithm. The distribution fo{R) is sampled after an initial transient, where the 
excited states components are filtered. All the sampling is accomplished iteratively through the integral equation 



/(i?,r)= J dR GdiR,R')GbiR,R')f{R',T-AT), 



(3.1) 



where 



Gd{R,R) = (47rDAT)-— exp 



{R-R -DArvoiR )f 
ADAt 



(3.2) 



vd = 2Vln*G, 



(3.3) 



Gb{R,R') = cxp{-{^)[El{R) + El{R')] + AtEt}, 



(3.4) 



here Et is a trial energy and E^ is the local energy given by 



f{R, t) for a long enough r goes to fo{R)- 



El = H^g/-^g] 



(3.5) 
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It is correct to write the Green's function as the product of Eq. 1)3. l|l only up to 0{t'^). This imphes that to obtain 
exact resuhs, within statistical fluctuations, short times At must be used in the iterations and an extrapolation to 
At —y performed (see however comments in Section llVII . This is the so called short-time approximation. 

Each configuration undergoes three steps: drift, diffusion and branching. Very frequently a single configuration 
is called a walker and an iteration of all walkers a generation. For the drift step we need to compute the quantum 
velocity vo. In the second step the configuration diffuses. This is accomplished by sampling Gd- Accordingly, a walker 
in R' is propagated during a time step At to its new point R through 



R = R' + X + DAtvd{R'), (3.6) 

where x ^-re normal deviates of a Gaussian function with variance 2_DAt and zero mean. 

To propagate a walker, we can change all the particle's coordinates at once or those of a single particle at a time. 
In this last case we perform N updates to propagate each walker. To improve^SiSS the approximation of the Green's 
function, we only accept moves with probability 



Pacccpt 



(i? R) 



1, 



Gd{R,R)^l{R) 

G<i(i?,i?')«'2j(i?') 



(3.7) 



and chose At such that more than 99% of the attempted moves are accepted. This condition imposes detailed balance 
on the splited Green's function and restores this property of the exact Green's function. Regardless of the time step, 
it guarantees also a correct sampling if hypothetically we could use as ipo the exact ground state wavefunction. In 
such case, this implementation of the DMC method reduces to a variational Monte Carlo calculation with trial moves 
sampled from Gd- 

To complete the iteration of Eq. we compute Gb considered as a weight for the walker. At the begin of 

the simulation, all the walker's weights are assumed to be equal to one. In order to minimize fluctuations in Gf,, an 
effective time step is usediSSiSl It is given by Ate// = At(Ap^/A/9^), where Ap^ is the mean square displacement of 
all proposed moves in the diffusion step and Ap^ is the related quantity when only accepted moves are considered. 
Finally the weight w' of the walker is updated to its new value w according to 



w = w'Gb{R,R'). (3.8) 

After propagating all walkers we have a new generation and a sample of the probability distribution / as an 
weighted average over the walkers. For reasons of efficiency the number of configurations used in the estimation of 
f{R) fluctuates according to the following branching rules. If w is greater than 2, the walker is duplicate and each 
one will carry half of its weight. On the other hand if two walkers, Ri and Rj, have weights less than 0.5, only one 
survives with a weight given by Wi + wj . The decision of which one will survive is made by sampling r = Wi/ {wi -\-Wj). 
That is, we draw a random number ^ and compare it with r. If r is small than ^ we keep configuration otherwise 
Ri is kept. If the weight computed by Eq. H3.8|l lies between 0.5 and 2, a single copy of the walker with weight w is 
kepted. The values of the weights where the branching is performed certainly could be changed, we only need rules 
that neither introduce bias or result in a scheme too inefficient. 

The number of walkers is controlled by adjusting the value of Et- This number is kept roughly constant. We have 
experienced with both heuristic and automatic changes of Et as given by 



Et = Eq + K\n{tp/cp), (3.9) 

where k is a parameter, tp is the target population and cp the current population. Adjustments in Et was made 
about once every 20 generations. For our purposes, the results were equivalent for both methods of changing Et- 



B. The Algorithm with Reweighting 



In many situations it is interesting to compute energies differences resulting from different interatomic potentials. 
However it is not always possible to simply use results from independent runs to obtain such differences. If they are 
small, statistical fluctuations might well produce errors that are bigger than these differences themselves, rendering 
the result meaningless. It is however possible to modify the DMC method in such a way that the same set of walkers 
are used to compute quantities of interest associated to the different potentials we want to investigate. The energies 
that are obtained are correlated and thus more meaningful differences can be computed. No approximations are 
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introduced. If the actual interest is on the energies, it is not necessary to use extrapolated estimators either. What 
we are proposing is to sample different probabilities distributions functions, associated to the different interatomic 
potentials we want to investigate, by using the same set of walkers with appropriate weights. As just mentioned, the 
values of the quantities of interest obtained are correlated and the errors associate with their difference reduced by 
orders of magnitude. We want to call attention to the fact that although our method relay in a set of weights it can 
not be related to the forward walking^S or rcptationpSS, methods. Always that we have a generation of equilibrated 
walkers, we can compute the quantities of interest without any further propagation. Moreover the weights we have 
for a given walker are associate to different interatomic potentials. 

In our modified DMC method, to each walker we attach a set of weights, one for each potential we want to consider. 
In our implementation of the algorithm we have attached three different weights for each walker, one for each of the 
three different interatomic potentials used. Of course this number in the method is arbitrary. It was chosen because 
of the specific aspects of many-body interactions in the interatomic potential we want to investigate. It would be 
possible to use only two weights or any other convenient number of weights. 

As in the standard algorithm, the calculations start with a set of walkers draw from l-fAcP- Since a single guide 
function is use, the drift and diffusion steps are performed exactly as before. We compute the drift velocity uu, 
generate the normal distribution of variates, update R according to Eq. (|3.()|) and accept it with probability paccept of 
Eq. 1)3. 7|l . In the present algorithm, we sample the three different probabilities distributions in which we are interested 
by completing the iteration considering three different Gb, and updating the weights as follows. 

To be specific let us consider a single walker just propagated to a new configuration R. One of its weights, it)'^' 

(2) 

is associate with the local energies E)^ (R) computed using only the Aziz two-body potential V2 in the Hamiltonian 
of Eq. H2.1|l . The weight w^^^ is updated according to Eq. (|3.8|l by evaluating Gb of Eq. (|3.4|) using the local energy 

Ef\R). Another weight of the same walker, w^'^^^ is calculated with the local energy E^^^\r) computed with 
the Hamiltonian that uses the V2D potential of Eq. 1)2. 2|l . it includes the triple-dipole contributions to the two- 
body potential. The calculation of the new value of w^"^^^ for this walker proceeds as before. Gb in Eq. (|3.4|) is 

evaluated employing E^^^\r) and the update finished according to Eq. H3.8|) . The third weight, w^^^"^), is associate 
to the interatomic potential that includes also exchange contributions. It is computed by considering E^^^^^{R) that 
depends on the Hamiltonian that uses the full potential V2DJ of Eq. (|2.6|1 . The update of w'-^^'-'^ is performed along 
the exact same lines already describe for the other weights. The different values of the weights are due only to the 
local energy used in their computation, i. e., to the interatomic potential employed. We remember again that the 
same configuration is used to compute these three weights. 

After propagating all walkers, a new generation will finally be obtained by using the following branching rules. When 
min(u>^^^ , u'^^^-' , ui^^^'^' ) is large than 2 the walker is duplicated and each one of the copies will carry half of the value 
of the weights: (u>(^^ /2, w'^-^^ /2, lo^^-^-^) /2). If two walkers, i and j have weights such that max(ti;j^^\ wf'^\ wf'^'^^) 
and max(u'j^', wj^^', wj^^'^') are less than 0.3, we consider each one of the weights individually. We draw a single 
random number ^ and make comparisons of ^ with r^'^\ j-C^D) ^^^^ ^(2DJ) ^ where r^^^ = wf'^ /{w^^'^ + '"^j''^)- Three 
situations might happen: i) in all the comparisons ^ is smaller than r'*"'', then wc keep walker i with weights {w^^^ + 
wf\ w^i^^^ +w^^^\ vji'^'^^ -l- u'j^^'^''} and discard walker j; ii) always ^ is greater than r^^\ in this case we keep Rj 
with the same sum of weights as above and discard Ri] iii) one of the comparisons favors a walker different from the 
other two. For definitencss let say that ^ is smaller than A"^^ and greater than A"^^^ and A'^^'^\ In this case we will 
keep the two walkers, Ri with weights {w-^^ -I- 0, 0} and Rj with weights {0, w'f^'' +w^^^\ w'f^'^'' -l-wj^'^''^}. 
These new weights are telling us that in fact we have deleted walker i from the calculations with the interatomic 
potential that includes three-body interactions and walker j when we are considering only the two-body potential. 
This is a bad situation in the sense that wc are introducing two walkers in the calculations that will not give anymore 
the correlations that we are looking for. Fortunately, if needed, the cases where this situation happens can be 
systematically reduced in a simple way. It is enough to decrease the threshold value used to combine walkers (see 
Section llV|) . Finally if one of the weights of a walker lies between 0.3 and 2, a single copy is maintained with weights 

It is useful to use only a single random number in the above comparisons. As expected each one of the calculations 
we are performing give exactly, within statistical fluctuations, the results obtained by the standard algorithm. Of 
course we could use three different random numbers in the comparison, however more frequently we would meet the 
unwanted situation described in iii) of the last paragraph. 

Periodically, about one every four or five generations wc compute several quantities of interest. Evaluations of the 

energies E^ ^ Em^^ and Em^'^'^ are readily obtained as weighted averages that include all the walkers Ri of the 
present generation: 

yZiW^^HRi) 
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Together with these quantities we have also evaluated the energy associate only to the damped triple-dipole term in 
the interatomic potential 

i?,^=i?r^-i?(^), (3.11) 

and the energy associate with the exchange term Vj 

Ei^E^^^^E^l (3.12) 
As already mentioned the computation of these values are straightforward because we have already estimates of the 

(k) 

exact energies Em , no extrapolations are needed. Along the runs, averages of these quantities arc formed and their 
estimates and associate errors obtained. 

IV. THE SIMULATIONS 

In the investigation of the properties of bulk helium we impose periodic boundary conditions. The cutoff convention, 
the distance beyond which a potential is set to zero, is enforced for all interactions at half of the box size, L/2. 
Distances between pairs of particles are computed by the minimum-image convention. When considering three-body 
interatomic interactions, the length of the third side of the triangles formed by the particles can not in general be 
computed using the minimum-image convention. A modification needs to be introduced so that the length of this 
side can be computed in a proper way and discarded if greater than L/2. To be specific let us consider particles «, j 
and k. We compute distances rij and rjfe using the minimum-image convention. The difference in the x coordinates 
of the associated particles are 

Xij — Xi Xj ^ij 1 

where the translation vector t is defined as 

Um = [{xi - Xm)/L]L 
and [x] is the closest integer to x. If the difference Xjk of the third side is computed as^i 

^jfc = X] - + Uj - Uk, (4.2) 

it is not hard to see that all possibilities in the simulation box arc taken into account and the right value of the 
distance can be obtained. If this value is little than L/2 the calculation for this triangle proceeds. For the y and z 
coordinates a similar approach is used and then the three-body interaction is computed if all sides for this triangle 
are lower than L/2. The calculation continues until all triangles have been considered. 

The diffusion Monte Carlo calculations started with an initial set of 400 configurations, previously draw from If/^Gp 
using the Metropolis algorithm. Before accumulating quantities of interest the excited states components of our 
ensemble of configurations are filtered by performing several iterations of Eq. I|3.1|l . This "equilibration" is typically 
of the order of 400 generations, and depends on the system density. 

We study the liquid phase using a guiding function of the Jastrow form 

i<j 

where the factor f{rij) = exp(— u(7'y )/2) explicitly correlates pairs of particles through a pseudopotential of the 
McMillan form u{rij) = (b/rij)^; 6 is a parameter. 
For the solid phase we have used a Nosanov-Jastrow 



guiding function, where 



^NjiR) = «'j(i?)$(i?) 



(4.4) 
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^R) = l[e^p\^ir.,^h)''\ (4.5) 

i 

is a mean field term that localizes the particles around the given lattice sites 1^. 

All guiding functions were previously optimized by performing variational calculations. Although this is a conve- 
nient way of obtaining the values of the parameters, in principle they could be obtained without pcrformating such 
calculations. It would be enough to chose parameters values that give the fastest filtering of the excited states of the 
initial configurations. 

The time steps At used in the calculations depends on the density. Their values vary within the range 0.001 to 
0.002 (K^^) in order to obtain more than 99% of acceptance of the attempted moves. We also observed that at this 
acceptance level, the extrapolation to At — > of the energies values were in excellent agreement, within statistical 
fluctuations, to the actual values obtained in the calculation itself. 

The quantities reported in this work were obtained by forming averages with about 500 estimates. Each estimate 
was performed after 4 generations. Blocking was used in order to avoid correlations in the calculations of the variances. 
The number of walkers during the simulations did not change by more than 10%. 

We have considered systems with 108 particles in the solid phase. In the liquid phase we have considered 64 particles. 
At po = 21.86 nm~^ to estimate size effects we have also performed simulations with 108 particles. Tail corrections of 
the two-body potential energy were made by assuming a pair distribution function equal to one beyond half the size of 
the simulation cell and integrating the potential up to infinity. No tail corrections were performed for the high-order 
interactions. For the Axilrod- Teller interaction, the tail correction is less than 7% of its value at po = 21.86 nm~'^ 
(see next section and Table II) . This value is in rough agreement with a previous estimate of this quantityi^ For the 
exchange energy the relative tail correction is bigger than the one of the dispersion energy. However it should remain 
within the statistical uncertainty of our results (see Tables I and II). 

The situation where we have a walkers with one of its weights equal to zero destroys the correlation we want to 
construct. If we combine walkers when all their weights is less or equal 0.3, we noticed that the number of walkers 
with at least one of weights equal zero does not exceed 2% of their total number. If needed this fraction can be 
further and systematically reduced by using a threshold smaller than 0.3 to combine walkers. As we have observed in 
our calculations, this is done at expense of a less efficient calculation. We have concluded that the threshold 0.3 for 
combination of walkers is perfectly reasonable for our purposes. 



V. RESULTS 



A. Liquid phase 

We conducted several independent runs at four different densities p of liquid helium, 19.64 nm~'^, at the experimental 
equilibrium density po = 21.86 nm~^, 24.01 nm~^ and at 26.23 nm~^. In Table I are shown the total energies obtained 
using the two-body potential V2 , the V2D potential of Eq. (|2.2|) , the V2 potential plus the Axilrod- Teller term, and the 
V2DJ potential of Ea. H2.6() . the exchange term added to the V2D potential. 

It is important to note that since the energies associated with these potentials arc calculated with a single set of 
walkers, they are correlated. We can believe that the results show their evolution as more elaborated interacting 
potentials are used, despite of the statistical uncertainties in the results. 

In Table II we shown very accurate calculations of the Axilrod- Teller and exchange energies at these four densities. 
For comparison we show also extrapolated results of perturbative calculations performed using configurations gen- 
erated with the V2 potential. We have plot these results in Figures 1 and 2. The Axilrod- Teller energies calculated 
using reweighting are greater than the extrapolated perturbative results. Moreover they do not always agree within 
the statistical uncertainty. The triplc-dipole interaction gives a positive contribution to the energy of the system and 
its value double when we go from the lowest to the highest density. 

The energies due to the exchange term in the V2DJ potential when calculated with respect to the total energies 
obtained with the V2D potential arc on average about 0.0010 K smaller than the extrapolated perturbative results. In 
addition there is no agreement within statistical uncertainties between the energies calculated with the reweighting 
and perturbative methods at p = 19.64 and p ~ 24.01 nm^'^. The energies calculated using reweighting are lower 
than the extrapolated perturbative results, contrary to what happens with the triple-dipole interaction. The exchange 
energy is also positive at all densities examined and increases with it. At the highest density it is approximately 50% 
greater than in the lowest one. 
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B. Solid phase 

For the solid phase we have considered four densities: 29.34, 32.88, 33.54 and 35.27 nm~'^. In Table I are shown the 
total energies, obtained with a systems of 108 particles in a fee structure. Again, we have considered the potentials 
1^2 1 and V2DJ- Table II shown our very accurate results of the Axilrod- Teller and exchange contributions to 
the potential energy and also extrapolated perturbative results for comparison. The difference is about the same we 
have observed in the liquid phase. For the Axilrod- Teller term they do not agree within the statistical uncertainty 
at p equal to 32.88 and 35.27 nm""^. In this phase as well, the Axilrod- Teller energies computed by reweighting are 
greater than the corresponding extrapolated perturbative results, and they remain positive. They also increase with 
the density. At the highest density (35.27 nm^^) it is 60% greater than in the lowest one. 

The contribution of the exchange term in the solid region is null or negative and differs significantly from the 
perturbative results. At the lowest density (29.34 nm~'^) the result obtained by reweighting gives a null contribution 
while the extrapolated perturbative quantity is positive. In the other densities the exchange energy computed by 
reweighting continues to be lower than the corresponding extrapolated perturbative values that are negative. In the 
solid phase the relative variation of the exchange energy is greater than the corresponding quantity for the damped 
Axilrod- Teller energy. 

C. Melting-Freezing Transition 

In order to follow the variations of the melting and freezing densities with the interacting potentials, we used a 
Maxwell (double-tangent) construction in analytical equations of state for the liquid and crystalline phases. The 
equations were determined by fits of our results to functions of the form 

E{p)^Eo + B{P^) + C{^^). (5.1) 
Po Po 

This functional form has been extensively used in the literature, including to fit experimental equation of stateiSiM 
We have fitted equations of state using results from the three different potentials, V2, V2D and ¥20,7- The fitted 
parameters, Eq, B, C and po in the liquid and solid phases are presented in Table III for these potentials. In Fig. 3 
we display results for the equations of state for the three potentials. 

The freezing and melting densities determined by the Maxwell double tangent construction are listed in Table IV. 
Looking at this table, we can follow the changes in the freezing and melting densities as more elaborated interacting 
potentials are used. The computed freezing densities differs, about 3% from the experimental value. This difference 
is of about 4% for the melting densities. The calculated freezing densities are below the experimental value, contrary 
to the computed melting densities that are above the experimental value. 

VI. DISCUSSION 

In this work we are able to verify without any approximations how small changes in the interacting potential affects 
some of the properties of bulk helium. It was possible to analyze how the Axilrod- Teller and the exchange three-body 
contributions to the interatomic potential modify the equation of state of this system. This was accomplished by 
using a single set of walkers with reweighting in a DMC calculation. The quantities of interest associate with the 
different potentials were obtained in a correlated fashion and so despite of the statistical errors their difference are 
meaningful. 

The total energies per atom presented in Table I show us that since our two-body potential is very accurate, high 
order terms in the description of the atomic interaction are needed. In the same token, as the contributions of the 
dispersion energy are much better known than those of the exchange energy, the results suggest that more efforts 
would be desirable in developing reliable ways of computing the energies associated to this last kind of interaction. 

The results of the Axilrod- Teller triple-dipole dispersion energy and of the three-body exchange energy as a function 
of the density reported in Table II, Figures 1 and 2, show qualitative agreement between the two methods we use in 
the calculations. However the usual approach of extrapolating perturbative calculations (made with configurations 
generated by the DMC method using a two-body potential) not always can be trusted in giving the right magnitude of 
this quantity. In some cases the results do not agree within statistical uncertainty with those obtained using reweight. 
This might be due to the simple functions used for extrapolation: functions of the Jastrow and Jastrow-Nosanov 
form, for the liquid and solid phases respectively. As expected, in both liquid and solid phases, the Axilrod- Teller 
term gives more important energy contribution to the total energy than the exchange term. Also it is interesting to 
note that at the lowest solid density we have considered (29.34 nm~'^), the extrapolated perturbative result associated 
with the exchange term gives a positive energy contribution whereas the value obtained by reweighting is null. 
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Although the mclusion of the three-body mteraction potentials used m this paper do not greatly modifies the melting 
and freezing densities, it is important to note that our calculations show that the inclusion of the triple-dipole and 
the exchange terms, as proposed by Cohen and Murell)^^ are leading the melting and freezing densities to their right 
values. All the small differences we observe in these quantities are within their estimate error (0.2 nm"'^). However 
since we used a single set of walkers to compute them, we can trust that their relative difference are significant in our 
calculations. Other consequence of including many body interactions in the potential can be seen in the calculation of 
the equilibrium density (Table III, Eq. (|OJ). For both potentials, V2D and V2DJ, the theoretical computed value of 
the equilibrium density becomes almost identical to the experimental value. The equilibrium density obtained using 
the full interatomic potential V2DJ diminished 0.3 nm~'^ from its value computed with the two-body potential. In the 
solid phase, in which the contribution of the exchange energy is greater, the parameter po decreases 0.75 nm^"^ in a 
similar comparison. 

The inclusion of three-body terms in the interatomic potential improves the agreement between computed and 
experimental values of properties like the binding energy, the equilibrium, melting and freezing densities. The results 
suggest that more reliable analytical expressions are needed for the calculation of the exchange energy in bulk helium. 
We could reach these conclusions by introducting reweighting in a DMC calculation, where a single set of walkers is 
used to compute properties associated to different potentials. Moreover we are able to calculate the energies associated 
to these potentials without extrapolation. 

The extension we have proposed to the DMC method might be very useful not only for the helium systems, but 
also for other quantum many-body systems where a clue is need to identify the best description between competing 
interacting potentials. Even if these potentials differ by very small amounts, the nature of the interactions they 
describe can be different. A better understanding of these differences and their relevance will enlarge our knowledge 
about the interactions and so about these systems themselves. The questions discussed in this work are not the only 
instance where the reweighting technique might be useful. The calculations of small energy contributions of spin-orbit 
terms in molecular physics^^ might be another situation where this method can help in a better understanding of a 
physical system. 
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TABLE I: Total energies per atom in units of K obtained at the given densities and potentials. Results in the second column 
for the potential V2, in the third and fourth columns the V2D and V2DJ potentials were considered (see text). In the liquid 
phase the results were obtained with 64 bodies and in the solid phase with 108 particles. In the last column we show the 
experimental values. 



p (nm ^) 




E 


(2) 
















D.J) 




Exp. 
















Liquid 














19.64 


-7 


121 


± 





006 


-7 


016 ±0 


006 


-7 


Oil 


±0 


006 


-7.01' 


21.86 


-7 


238 


± 





009 


-7 


103 ±0 


010 


-7 


097 


±0 


010 


-7.14** 


21.86" 


-7 


240 


± 





007 


-7 


101 ±0 


007 


-7 


095 


±0 


007 


-7.14' 


24.01 


-7 


120 


± 





010 


-6 


949 ±0 


010 


-6 


942 


±0 


010 


-7.00' 


26.23 


-6 


541 


± 





014 


-6 


325 ±0 
Solid 


014 


-6 


318 


±0 


014 


-6.53' 


29.34 


-5 


907 


± 





004 


-5 


600 ±0 


003 


-5 


600 


±0 


004 


-5.78" 


32.88 


-4 


489 


± 





006 


-4 


071 ±0 


006 


-4 


076 


±0 


007 




33.54 


-4 


089 


± 





005 


-3 


648 ±0 


005 


-3 


656 


±0 


005 


-3.94" 


35.27 


-2 


831 


± 





006 


-2 


323 ±0 


006 


-2 


336 


±0 


006 


-2.70" 



"Result for 108 particles. 
'Reference \ 3A 
"Reference l3fiL 
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TABLE II: Energies per particle in units of K associated to the triple-dipole term (E^) and the exchange term {E'') at the 
given densities obtained by reweighting and by extrapolation of perturbative calculation. 







E 


7 


P (nm ^) 


Rew. 


Extr. 


Rew. 


Extr. 






Liquid 






19.64 


0.105 ±0.001 


0.1012 ±0.0002 


0.0044 ± 0.0003 


0.0056 ± 0.0001 


21.86 


0.135 ±0.001 


0.1333 ±0.0002 


0.0056 ± 0.0004 


0.0066 ± 0.0001 


21.86" 


0.139 ±0.002 


0.1351 ±0.0006 


0.0058 ± 0.0004 


0.0068 ± 0.0001 


24.01 


0.171 ±0.001 


0.1698 ± 0.0002 


0.0063 ± 0.0003 


0.0074 ± 0.0001 


26.23 


0.217 ±0.001 


0.2136 ±0.0002 


0.0069 ± 0.0005 


0.0078 ± 0.0001 






Sohd 






29.34 


0.307 ±0.001 


0.3035 ± 0.0002 


0.0000 ± 0.0003 


0.0019 ± 0.0001 


32.88 


0.418 ±0.001 


0.4140 ± 0.0002 


-0.0051 ± 0.0007 


-0.0028 ± 0.0002 


33.54 


0.441 ±0.001 


0.4379 ± 0.0002 


-0.0072 ± 0.0008 


-0.0050 ± 0.0002 


35.27 


0.508 ±0.001 


0.5029 ± 0.0002 


-0.0128 ± 0.0006 


-0.0097 ± 0.0002 



"Result for 108 particles. 
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TABLE III: Fitting parameters of the liquid and solid equations of state for three different potentials. In the first line, for 
both the liquid and solid phase, the two-body potential V2 of Aziz et al. (Ref. 0) was used. Then we present results when the 
three-body Axilrod- Teller term is included in the interacting potential, V2D. In the rows with V2DJ we show results obtained 
when the full potential, that includes the three-body exchange term in V2D, was used. The units of -Bo, B and C are expressed 
in K. 



Potential po (nm 




B 


C 






Liquid 






V2 


22.133 


-7.240 


13.549 


37.025 


V2D 


21.845 


-7.103 


12.143 


35.705 


V2DJ 


21.834 


-7.097 


12.081 


35.480 






Solid 






V2 


26.795 


-6.200 


31.880 


5.661 


V2D 


26.399 


-5.980 


29.739 


7.870 


V2D.J 


26.045 


-6.028 


25.233 


11.844 
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TABLE IV: Melting and freezing densities using three different potentials calculated by the Maxwell double tangent construction 
method. The V2 potential is the two-body potential of Aziz et al. (Ref. 0)- The V2D potential is build using the V2 potential 
plus the three-body Axilrod- Teller interaction term. Finally the V2DJ potential includes the V2 potential, the three-body 
Axilrod- Teller and the exchange terms. In the last line we give the experimental values. 



Potential 


Pi (mn^^) 


Pm (nm ^) 


V2 


24.94 


29.39 


V2D 


25.00 


29.35 


V2D.J 


24.99 


29.28 


Exp. 


25.8" 


28.0" 



"Reference Isvi 
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CAPTIONS 

Figure 1.- Energy per atom associated to the three-body AxUrod- Teller interaction term for the liquid and solid 
phases. The crosses stand for the DMC results with rewcighting. The circles show extrapolated estimates of pertur- 
bative calculations. The statistical errors are smaller than the size of the symbols. The results were obtained using a 
simulation cell with 64 particles for the liquid phase and 108 for the solid one. 

Figure 2.- Energy per particle associated to the three-body exchange term for the liquid and solid phases, DMC 
with rewcighting (crosses) and extrapolated perturbative results (circles) . The statistical errors of the last calculations 
are smaller than the size of the symbols. The results were obtained using a simulation cell with 64 particles for the 
liquid phase and 108 for the solid one. 

Figure 3.- Analytical equations of state with three different potentials for the solid and liquid phases. The doted 
line represent the equation of state obtained using the results determined with the two-body potential V2. The solid 
line represent results using the V2D and V2DJ potentials that includes only the Axilrod- Teller term and this term plus 
the exchange one, respectively. At the figure scale the two fits are indistinguishable. The points represent results 
from our calculations. 




Density (nm" ) 




Density (nm ) 




Density (nm ) 



